Based on the alpha and beta diversity results, I am not expecting to find a huge difference between groups in terms of bacterial composition. Please note that many of these code are from the Microbiome Intensive Course (MIC Course, Duke University, 2023) and thus the work of Dr. Pixu Shi.

This is after talking with Dr. Shi about what analyses I want to run.

1 Set-Up

First, update the ANCOMBC package if necessary.

#if (!requireNamespace("BiocManager", quietly=TRUE))
#     install.packages("BiocManager")
#BiocManager::install("ANCOMBC")

2 Load libraries and directories

# rm(list = ls(all = TRUE)) # Unload all packages if necessary
library(phyloseq)
library(fs)
library(rstatix)
library(ggpubr)
library(microViz)
library(ANCOMBC)
library(pals)
library(RColorBrewer)
library(colorBlindness)
library(microbiome)
library(ggrepel)
library(DT)
library(plotly)
library(DESeq2)
library(parallel)
library(dplyr)
library(tidyr)
library(readr)
library(tibble)
library(stringr)
rm(list=ls()) # clear environment 

git.dir = file.path("/hpc/group/dmc/Eva/GLP_KO_Microbiome_final") # CHANGE ME 
ps.rds <- file.path(git.dir, "ps.GLP1.rds")
ps <- read_rds(ps.rds)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 234 samples ]
## sample_data() Sample Data:       [ 234 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq()      DNAStringSet:      [ 3304 reference sequences ]
fig.dir = file.path(git.dir, "figures")

map.file = file.path(git.dir, "Metadata_GLP1.csv")
meta.df = read_csv(map.file, show_col_types = FALSE)

3 Prepare phyloseq object

I will now aggregate at the genus level, as it is the smallest taxonomic rank that analyses will be run on. In other words, I do not plan on running differential abundance analysis at ASV level.

Let’s prepare the phyloseq object for differential abundance analysis (DAA). Since we are interested in 1) differentially abundant taxa and 2) correlation analysis, we will limit the taxa to only ones that appear in > 3 samples for each sample type.

To avoid running analyses on extremely rare taxa that may be sequencing artifacts, the threshold for minimum number of counts is set at 10. This is, of course, arbitrary. 10 was chosen as 1% of number of read threshold of 1000 above for pruning samples. Based on how I wrote the code, taxa must be present have > 10 reads or more in > 3 samples.

Furthermore, categorical variables should be re-leveled as factors.

# Subset to true fecal samples only and remove tree
ps  %>% subset_samples(Sample_type == "Fecal" & Type == "True Sample") -> ps_fecal
ps_fecal_notree <- phyloseq(otu_table(ps_fecal), tax_table(ps_fecal), sample_data(ps_fecal))

# Categorical variables should be factors 
ps_fecal_notree@sam_data$Genotype <- factor(ps_fecal_notree@sam_data$Genotype, levels = c("WT", "KO"))
ps_fecal_notree@sam_data$Sex <- factor(ps_fecal_notree@sam_data$Sex, levels = c("Female", "Male"))
ps_fecal_notree@sam_data$Intranasal_Treatment <- factor(ps_fecal_notree@sam_data$Intranasal_Treatment, levels = c("PBS", "HDM"))
ps_fecal_notree@sam_data$Surgery <- factor(ps_fecal_notree@sam_data$Surgery, levels = c("None", "Sham", "VSG"))
ps_fecal_notree@sam_data$Mouse %>% unique() -> mouse_levels
ps_fecal_notree@sam_data$Mouse <- factor(ps_fecal_notree@sam_data$Mouse, levels = mouse_levels)
ps_fecal_notree@sam_data$Diet <- factor(ps_fecal_notree@sam_data$Diet, levels = c("Normal_Chow", "High_Fat_Diet"))
ps_fecal_notree@sam_data$Group <- factor(ps_fecal_notree@sam_data$Group, levels = c("NA", "Control", "1", "2"))
ps_fecal_notree@sam_data$Week <- factor(ps_fecal_notree@sam_data$Week, levels = c("0", "10", "13"))
ps_fecal_notree
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3304 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 3304 taxa by 6 taxonomic ranks ]
# Aggregate at genus level
ps_g_notree_fecal <- tax_glom(ps_fecal_notree, "Genus")
taxa_names(ps_g_notree_fecal) <- tax_table(ps_g_notree_fecal)[, 'Genus']

ps_g_notree_fecal
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 314 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 314 taxa by 6 taxonomic ranks ]

4 Getting ready: prune data

We will prune the phyloseq object at each phylogenetic rank at and above genus. Since we plan on running correlation analysis after differential abundance, let’s test for taxa that are present in at least 5 samples (arbitrary number) with greater than 10 reads in each sample (also arbitrary. It’s also 1% of 1000 reads, which was the threshold used to remove samples from analysis).

x` Genus

# Pruning: set threshold number 
threshold_num = 4
nreads_prune = 10

# Prune taxa at each phylogenetic rank 
ps_g_notree_fecal %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_genus.prune
ps_genus.prune
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 70 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 70 taxa by 6 taxonomic ranks ]
ps_g_notree_fecal %>% transform_sample_counts(function(x) x/sum(x)) -> ps_genus.ts
ps_genus.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 314 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 314 taxa by 6 taxonomic ranks ]
subset_taxa(ps_genus.ts, ps_genus.ts@tax_table[, "Genus"] %in% ps_genus.prune@tax_table[, "Genus"]) -> ps_genus.prune.ts
ps_genus.prune.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 70 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 70 taxa by 6 taxonomic ranks ]
rowSums(ps_genus.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> genus_prune_RelAb

summary(genus_prune_RelAb$RelAb)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   93.13   99.92   99.97   99.89  100.00  100.00

4.1 Family

ps_fecal_notree %>% tax_glom("Family") -> ps_Family  
taxa_names(ps_Family) <- tax_table(ps_Family)[, 'Family']
ps_Family
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 169 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 169 taxa by 6 taxonomic ranks ]
ps_Family %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Family.ts
ps_Family.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 169 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 169 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank 
ps_Family %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Family.prune
ps_Family.prune
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 37 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 37 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Family.ts, ps_Family.ts@tax_table[, "Family"] %in% ps_Family.prune@tax_table[, "Family"]) -> ps_Family.prune.ts
ps_Family.prune.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 37 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 37 taxa by 6 taxonomic ranks ]
rowSums(ps_Family.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Family_prune_RelAb

summary(Family_prune_RelAb)
##      RelAb       
##  Min.   : 94.21  
##  1st Qu.: 99.99  
##  Median :100.00  
##  Mean   : 99.95  
##  3rd Qu.:100.00  
##  Max.   :100.00

4.2 Order

ps_fecal_notree %>% tax_glom("Order") -> ps_Order  
taxa_names(ps_Order) <- tax_table(ps_Order)[, 'Order']
ps_Order
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 109 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 109 taxa by 6 taxonomic ranks ]
ps_Order %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Order.ts
ps_Order.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 109 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 109 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank 
ps_Order %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Order.prune
ps_Order.prune
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 26 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 26 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Order.ts, ps_Order.ts@tax_table[, "Order"] %in% ps_Order.prune@tax_table[, "Order"]) -> ps_Order.prune.ts
ps_Order.prune.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 26 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 26 taxa by 6 taxonomic ranks ]
rowSums(ps_Order.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Order_prune_RelAb

summary(Order_prune_RelAb)
##      RelAb       
##  Min.   : 94.22  
##  1st Qu.:100.00  
##  Median :100.00  
##  Mean   : 99.95  
##  3rd Qu.:100.00  
##  Max.   :100.00

4.3 Class

ps_fecal_notree %>% tax_glom("Class") -> ps_Class  
taxa_names(ps_Class) <- tax_table(ps_Class)[, 'Class']
ps_Class
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 42 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 42 taxa by 6 taxonomic ranks ]
ps_Class %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Class.ts
ps_Class.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 42 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 42 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank 
ps_Class %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Class.prune
ps_Class.prune
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 12 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 12 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Class.ts, ps_Class.ts@tax_table[, "Class"] %in% ps_Class.prune@tax_table[, "Class"]) -> ps_Class.prune.ts
ps_Class.prune.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 12 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 12 taxa by 6 taxonomic ranks ]
rowSums(ps_Class.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Class_prune_RelAb

summary(Class_prune_RelAb)
##      RelAb       
##  Min.   : 99.99  
##  1st Qu.:100.00  
##  Median :100.00  
##  Mean   :100.00  
##  3rd Qu.:100.00  
##  Max.   :100.00

4.4 Phylum

ps_fecal_notree %>% tax_glom("Phylum") -> ps_Phylum  
taxa_names(ps_Phylum) <- tax_table(ps_Phylum)[, 'Phylum']
ps_Phylum
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]
ps_Phylum %>% transform_sample_counts(function(x) x/sum(x)) -> ps_Phylum.ts
ps_Phylum.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]
# Prune taxa at each phylogenetic rank 
ps_Phylum %>%
  filter_taxa(., function(x) {sum(x > nreads_prune) > threshold_num}, prune = TRUE) -> ps_Phylum.prune
ps_Phylum.prune
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 6 taxonomic ranks ]
subset_taxa(ps_Phylum.ts, ps_Phylum.ts@tax_table[, "Phylum"] %in% ps_Phylum.prune@tax_table[, "Phylum"]) -> ps_Phylum.prune.ts
ps_Phylum.prune.ts
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 178 samples ]
## sample_data() Sample Data:       [ 178 samples by 42 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 6 taxonomic ranks ]
rowSums(ps_Phylum.prune.ts@otu_table) %>% as.data.frame() %>% 
  dplyr::rename(RelAb = ".") %>% mutate(RelAb = RelAb * 100,
                                        RelAb = format(RelAb, scientific = FALSE),
                                        RelAb = as.numeric(RelAb)) %>% 
  dplyr::arrange(RelAb) -> Phylum_prune_RelAb

summary(Phylum_prune_RelAb)
##      RelAb    
##  Min.   :100  
##  1st Qu.:100  
##  Median :100  
##  Mean   :100  
##  3rd Qu.:100  
##  Max.   :100

5 Visualization

Let’s visualize at Phylum level to see how much information was kept:

plot_bar(ps_genus.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Genus", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Family.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Family", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Order.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Order", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Class.prune.ts, x = "Label", fill="Class") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Class, fill=Class), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Class", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

plot_bar(ps_Phylum.prune.ts, x = "Label", fill="Phylum") + 
  facet_grid(scales="free", space = "free_x") + 
  geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +  
  theme_bw()  +   
  coord_cartesian(ylim = c(0,1), expand=0) +
  theme(axis.text.x=element_blank()) + # not showing samples names 
  labs(y="Relative Abundance", title = "Relative Abundance of fecal Samples, Pruned @Phylum", x = "Sample")+
  scale_color_manual(values=SteppedSequential5Steps) + 
  scale_fill_manual(values=SteppedSequential5Steps) 

6 Subset metadata

meta.df %>% filter(Type == "True Sample" & Sample_type == "Fecal") %>% 
  filter(Label %in% phyloseq::sample_data(ps_g_notree_fecal)$Label) -> meta.fecal
meta.fecal

7 Q1: Fecal microbiome, HFD, wk 10 & 13, Surgery * Week

Let’s see how many mice we have:

meta.fecal %>% filter(Diet == "High_Fat_Diet" & Week != 0) -> meta.f.hfd
meta.f.hfd %>% 
  dplyr::count(Week, Surgery, Genotype) %>% 
  arrange(n)
meta.f.hfd %>% 
  dplyr::count(Week, Surgery) %>% 
  arrange(n)

I will run Week * Surgery.

nrow(meta.f.hfd)
## [1] 78

7.1 Make TSE objects

ps_genus.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Genus
tse.fecal.Q1.Genus
## class: TreeSummarizedExperiment 
## dim: 70 78 
## metadata(0):
## assays(1): counts
## rownames(70): Peptococcus Anaeroplasma ... Bifidobacterium
##   Corynebacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Family.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Family
tse.fecal.Q1.Family
## class: TreeSummarizedExperiment 
## dim: 37 78 
## metadata(0):
## assays(1): counts
## rownames(37): Acholeplasmataceae Peptococcaceae ... Bifidobacteriaceae
##   Corynebacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Order.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Order
tse.fecal.Q1.Order
## class: TreeSummarizedExperiment 
## dim: 26 78 
## metadata(0):
## assays(1): counts
## rownames(26): Acholeplasmatales Peptococcales ... Bifidobacteriales
##   Corynebacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Class.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Class
tse.fecal.Q1.Class
## class: TreeSummarizedExperiment 
## dim: 12 78 
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
##   Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Phylum.prune %>% subset_samples(Diet == "High_Fat_Diet" & Week !="0")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q1.Phylum
tse.fecal.Q1.Phylum
## class: TreeSummarizedExperiment 
## dim: 10 78 
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
##   Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(78): SRR29054904 SRR29054909 ... SRR29055192 SRR29055193
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

7.2 ANCOM2

7.2.1 Genus

set.seed(123) # Set seed for reproducibility 
n_cl_requested = 24

# Genus level
fecal.Q1.Genus <- ANCOMBC::ancom(data=tse.fecal.Q1.Genus, assay_name="counts", tax_level="Genus",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Week",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Genus = fecal.Q1.Genus$res
res.fecal.Q1.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Genus.sig
tab.fecal.Q1.Genus.sig 

7.3 Family

set.seed(123) # Set seed for reproducibility 

# Family level
fecal.Q1.Family <- ANCOMBC::ancom(data=tse.fecal.Q1.Family, assay_name="counts", tax_level="Family",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Week",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Family = fecal.Q1.Family$res
res.fecal.Q1.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Family.sig
tab.fecal.Q1.Family.sig 

7.3.1 Order

set.seed(123) # Set seed for reproducibility 

# Order level
fecal.Q1.Order <- ANCOMBC::ancom(data=tse.fecal.Q1.Order, assay_name="counts", tax_level="Order",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Week",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Order = fecal.Q1.Order$res
res.fecal.Q1.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Order.sig
tab.fecal.Q1.Order.sig 

7.3.2 Class

set.seed(123) # Set seed for reproducibility 

# Class level
fecal.Q1.Class <- ANCOMBC::ancom(data=tse.fecal.Q1.Class, assay_name="counts", tax_level="Class",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Week",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Class = fecal.Q1.Class$res
res.fecal.Q1.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Class.sig
tab.fecal.Q1.Class.sig 

7.3.3 Phylum

set.seed(123) # Set seed for reproducibility 

# Phylum level
fecal.Q1.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q1.Phylum, assay_name="counts", tax_level="Phylum",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Surgery", adj_formula="Week",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q1.Phylum = fecal.Q1.Phylum$res
res.fecal.Q1.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q1.Phylum.sig
tab.fecal.Q1.Phylum.sig 

7.3.4 Aggregate Results

tab.fecal.Q1.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q1.Family.sig$Tax_rank <- "Family"
tab.fecal.Q1.Order.sig$Tax_rank <- "Order"
tab.fecal.Q1.Class.sig$Tax_rank <- "Class"
tab.fecal.Q1.Phylum.sig$Tax_rank <- "Phylum"

rbind(tab.fecal.Q1.Genus.sig,
      tab.fecal.Q1.Family.sig,
      tab.fecal.Q1.Order.sig,
      tab.fecal.Q1.Class.sig,
      tab.fecal.Q1.Phylum.sig) -> tab.fecal.Q1

tab.fecal.Q1$Test <- "ANCOM II, fecal microbiome: all HFD, wk 10 & 13, Surgery * Week. Structural zeroes from Surgery"

tab.fecal.Q1 

7.4 DESeq2

7.4.1 Genus

ds.fecal.Q1_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q1.Genus, ~ Surgery * Week)
## converting counts to integer mode
dds.fecal.Q1_Genus <- DESeq2::DESeq(ds.fecal.Q1_Genus)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 8 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Genus)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1015262       2530      52223        12554.      13016.        6735.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    70          0           39.2
resultsNames(dds.fecal.Q1_Genus)
## [1] "Intercept"            "Surgery_Sham_vs_None" "Surgery_VSG_vs_None" 
## [4] "Week_13_vs_10"        "SurgerySham.Week13"   "SurgeryVSG.Week13"
results(dds.fecal.Q1_Genus, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Genus.Res
fecal.Q1.Genus.Res$Comparison <- c("Week: week 13 over week 10")

results(dds.fecal.Q1_Genus, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res

results(dds.fecal.Q1_Genus, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res

results(dds.fecal.Q1_Genus, name  = c("SurgerySham.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery Sham, wk 13")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res

results(dds.fecal.Q1_Genus, name  = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Genus.Res) -> fecal.Q1.Genus.Res

fecal.Q1.Genus.Res$Compare_taxon <- "Genus"

fecal.Q1.Genus.Res

7.4.2 Family

ds.fecal.Q1_Family <- DESeq2::DESeqDataSet(tse.fecal.Q1.Family, ~ Surgery * Week)
## converting counts to integer mode
dds.fecal.Q1_Family <- DESeq2::DESeq(ds.fecal.Q1_Family)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Family)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1280900       3584      68847        15486.      16422.        8868.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    37          0           22.8
resultsNames(dds.fecal.Q1_Family)
## [1] "Intercept"            "Surgery_Sham_vs_None" "Surgery_VSG_vs_None" 
## [4] "Week_13_vs_10"        "SurgerySham.Week13"   "SurgeryVSG.Week13"
results(dds.fecal.Q1_Family, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Family.Res
fecal.Q1.Family.Res$Comparison <- c("Week: week 13 over week 10")

results(dds.fecal.Q1_Family, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Family.Res) -> fecal.Q1.Family.Res

results(dds.fecal.Q1_Family, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Family.Res) -> fecal.Q1.Family.Res

results(dds.fecal.Q1_Family, name  = c("SurgerySham.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery Sham, wk 13")
rbind(tmp, fecal.Q1.Family.Res) -> fecal.Q1.Family.Res

results(dds.fecal.Q1_Family, name  = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) 
fecal.Q1.Family.Res$Compare_taxon <- "Family"

fecal.Q1.Family.Res

7.4.3 Order

ds.fecal.Q1_Order <- DESeq2::DESeqDataSet(tse.fecal.Q1.Order, ~ Surgery * Week)
## converting counts to integer mode
dds.fecal.Q1_Order <- DESeq2::DESeq(ds.fecal.Q1_Order)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 4 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Order)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1286962       3587      68977        15488.      16500.        8913.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    26          0           15.2
resultsNames(dds.fecal.Q1_Order)
## [1] "Intercept"            "Surgery_Sham_vs_None" "Surgery_VSG_vs_None" 
## [4] "Week_13_vs_10"        "SurgerySham.Week13"   "SurgeryVSG.Week13"
results(dds.fecal.Q1_Order, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Order.Res
fecal.Q1.Order.Res$Comparison <- c("Week: week 13 over week 10")

results(dds.fecal.Q1_Order, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res

results(dds.fecal.Q1_Order, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res

results(dds.fecal.Q1_Order, name  = c("SurgerySham.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery Sham, wk 13")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res

results(dds.fecal.Q1_Order, name  = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Order.Res) -> fecal.Q1.Order.Res

fecal.Q1.Order.Res$Compare_taxon <- "Order"

fecal.Q1.Order.Res

7.4.4 Class

ds.fecal.Q1_Class <- DESeq2::DESeqDataSet(tse.fecal.Q1.Class, ~ Surgery * Week)
## converting counts to integer mode
dds.fecal.Q1_Class <- DESeq2::DESeq(ds.fecal.Q1_Class)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Class)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1288152       3587      69016        15488.      16515.        8912.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    12          0           9.33
resultsNames(dds.fecal.Q1_Class)
## [1] "Intercept"            "Surgery_Sham_vs_None" "Surgery_VSG_vs_None" 
## [4] "Week_13_vs_10"        "SurgerySham.Week13"   "SurgeryVSG.Week13"
results(dds.fecal.Q1_Class, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Class.Res
fecal.Q1.Class.Res$Comparison <- c("Week: week 13 over week 10")

results(dds.fecal.Q1_Class, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: Sham over None")
rbind(tmp, fecal.Q1.Class.Res) -> fecal.Q1.Class.Res

results(dds.fecal.Q1_Class, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Class.Res) -> fecal.Q1.Class.Res

results(dds.fecal.Q1_Class, name  = c("SurgerySham.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) 
results(dds.fecal.Q1_Class, name  = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Class.Res) -> fecal.Q1.Class.Res

fecal.Q1.Class.Res$Compare_taxon <- "Class"

fecal.Q1.Class.Res

7.4.5 Phylum

ds.fecal.Q1_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q1.Phylum, ~ Surgery * Week)
## converting counts to integer mode
dds.fecal.Q1_Phylum <- DESeq2::DESeq(ds.fecal.Q1_Phylum)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 1 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
DESeq2::summary(dds.fecal.Q1_Phylum)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1288193       3587      69016        15490.      16515.        8912.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    10          0              8
resultsNames(dds.fecal.Q1_Phylum)
## [1] "Intercept"            "Surgery_Sham_vs_None" "Surgery_VSG_vs_None" 
## [4] "Week_13_vs_10"        "SurgerySham.Week13"   "SurgeryVSG.Week13"
results(dds.fecal.Q1_Phylum, tidy = TRUE, contrast = c("Week", "13", "10")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q1.Phylum.Res
fecal.Q1.Phylum.Res$Comparison <- c("Week: week 13 over week 10")

results(dds.fecal.Q1_Phylum, contrast = c("Surgery", "Sham", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) 
results(dds.fecal.Q1_Phylum, contrast = c("Surgery", "VSG", "None"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery: VSG over None")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res

results(dds.fecal.Q1_Phylum, name  = c("SurgerySham.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) 
results(dds.fecal.Q1_Phylum, name  = c("SurgeryVSG.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Surgery VSG, wk 13")
rbind(tmp, fecal.Q1.Phylum.Res) -> fecal.Q1.Phylum.Res

fecal.Q1.Phylum.Res$Compare_taxon <- "Phylum"

fecal.Q1.Phylum.Res

7.4.6 Aggregate results

rbind(fecal.Q1.Genus.Res, 
      fecal.Q1.Family.Res,
      fecal.Q1.Class.Res,
      fecal.Q1.Order.Res,
      fecal.Q1.Phylum.Res) -> fecal.Q1.deseq

fecal.Q1.deseq$Test <- "DESeq2, Fecal microbiome: all HFD, wk 10 + 13, Surgery * Week"

fecal.Q1.deseq

7.5 Intersection

tab.fecal.Q1 %>% filter(taxon %in% fecal.Q1.deseq$row) -> ancom.intersect.hfd
ancom.intersect.hfd
fecal.Q1.deseq %>% filter(row %in% tab.fecal.Q1$taxon) 

7.6 Graph

psmelt(ps_genus.prune.ts) -> Genus.prune.melt 
psmelt(ps_Family.prune.ts) -> Family.prune.melt
psmelt(ps_Order.prune.ts) -> Order.prune.melt
psmelt(ps_Class.prune.ts) -> Class.prune.melt
psmelt(ps_Phylum.prune.ts) -> Phylum.prune.melt
bind_rows(
  Genus.prune.melt,
  Family.prune.melt,
  Order.prune.melt,
  Class.prune.melt,
  Phylum.prune.melt
)  -> fecal.RA.Full

fecal.RA.Full$Surgery <- factor(fecal.RA.Full$Surgery, levels = c("None", "Sham", "VSG"))

fecal.RA.Full %>% 
  dplyr::select(OTU:Label, Sample_type:Mouse, Genotype:Surgery, Week) -> fecal.RA.Full

head(fecal.RA.Full)
tail(fecal.RA.Full)
fecal.RA.Full %>% filter(Diet == "High_Fat_Diet" & Week != "0") -> fecal.graph.hfd
head(fecal.graph.hfd)
tail(fecal.graph.hfd)
surgery.labs <- c("No Surgery", "Sham Surgery", "VSG")
names(surgery.labs) <- c("None", "Sham", "VSG")

fecal.graph.hfd %>% 
  filter(OTU %in% ancom.intersect.hfd$taxon) %>%   
  ggplot(aes(x = Surgery, y = Abundance, color = Surgery)) + 
  geom_boxplot() + geom_point() +  
  labs(y="Relative Abundance", title = str_glue(ancom.intersect.hfd$taxon[1], " Relative Abundance, fecal microbiome"),
       x="Surgery") + 
  theme_bw(base_size = 14) + 
  facet_wrap(~Week)

fecal.graph.hfd %>% 
  filter(OTU %in% ancom.intersect.hfd$taxon) %>%   
  ggplot(aes(x = Week, y = Abundance, color = Week)) + 
  geom_boxplot() + geom_point() +  
  labs(y="Relative Abundance", title = str_glue(ancom.intersect.hfd$taxon[1], " Relative Abundance, fecal microbiome"),
       x="Surgery") + 
  theme_bw(base_size = 14) + 
  facet_wrap(~Surgery, scales = "free")

fecal.graph.hfd %>% 
  filter(OTU %in% ancom.intersect.hfd$taxon) %>% 
  group_by(Week) %>%
  kruskal_test(Abundance ~ Surgery) %>%
  adjust_pvalue(method = "holm") %>%
  add_significance()

Post-hoc: pairwise comparison with Wilcoxon test.

fecal.graph.hfd %>% 
  filter(OTU %in% ancom.intersect.hfd$taxon) %>% 
  group_by(Week) %>% 
  wilcox_test(Abundance ~ Surgery, p.adjust.method = "holm") %>%
  add_significance() %>% 
  add_y_position() -> fecal.hfd.res

fecal.hfd.res
week.labs <- c("Week 0", "Week 10", "Week 13")
names(week.labs) <- c("0", "10", "13")

fecal.graph.hfd %>% 
  filter(OTU %in% ancom.intersect.hfd$taxon) %>%   
  ggplot(aes(x = Surgery, y = Abundance, color = Surgery)) + 
  geom_boxplot() + geom_point() +  
  labs(y="Lactobacillaceae Relative Abundance",
       x="Surgery") + 
  theme_bw(base_size = 14) + 
  facet_wrap(~Week,  labeller = labeller(Week = week.labs)) + theme(legend.position = "top") -> fecal.hfd.res.graph

fecal.hfd.res.graph

fecal.hfd.res$y.position <- c(0.45, 0.5, 0.57, 0.55, 0.64, 0.7)
fecal.hfd.res.graph  + 
  stat_pvalue_manual(fecal.hfd.res, label = "p.adj.signif") -> fecal.hfd.res.graph

fecal.hfd.res.graph

ragg::agg_jpeg(file.path(fig.dir, "fecal_surgery_lactobacillaceae.jpeg"), width = 7, height = 5, units = "in", res = 600)
fecal.hfd.res.graph
dev.off()
## png 
##   2

8 Q2: no surgery group, all weeks, Diet * Week + Genotype

Let’s see how many mice we have:

meta.fecal %>% filter(Surgery == "None") -> meta.f.nosurgery
meta.f.nosurgery %>% 
  dplyr::count(Week, Diet, Genotype) %>% 
  arrange(n)
meta.f.nosurgery %>% 
  dplyr::count(Week, Diet) %>% 
  arrange(n)

I will run Diet * Week + Genotype.

nrow(meta.f.nosurgery)
## [1] 120

8.1 Make TSE objects

ps_genus.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Genus
tse.fecal.Q2.Genus
## class: TreeSummarizedExperiment 
## dim: 70 120 
## metadata(0):
## assays(1): counts
## rownames(70): Peptococcus Anaeroplasma ... Bifidobacterium
##   Corynebacterium
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Family.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Family
tse.fecal.Q2.Family
## class: TreeSummarizedExperiment 
## dim: 37 120 
## metadata(0):
## assays(1): counts
## rownames(37): Acholeplasmataceae Peptococcaceae ... Bifidobacteriaceae
##   Corynebacteriaceae
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Order.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Order
tse.fecal.Q2.Order
## class: TreeSummarizedExperiment 
## dim: 26 120 
## metadata(0):
## assays(1): counts
## rownames(26): Acholeplasmatales Peptococcales ... Bifidobacteriales
##   Corynebacteriales
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Class.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Class
tse.fecal.Q2.Class
## class: TreeSummarizedExperiment 
## dim: 12 120 
## metadata(0):
## assays(1): counts
## rownames(12): Campylobacteria Desulfovibrionia ... Clostridia
##   Actinobacteria
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL
ps_Phylum.prune %>% subset_samples(Surgery == "None")  %>%
  mia::makeTreeSummarizedExperimentFromPhyloseq() -> tse.fecal.Q2.Phylum
tse.fecal.Q2.Phylum
## class: TreeSummarizedExperiment 
## dim: 10 120 
## metadata(0):
## assays(1): counts
## rownames(10): Campylobacterota Desulfobacterota ... Firmicutes
##   Actinobacteriota
## rowData names(6): Kingdom Phylum ... Family Genus
## colnames(120): SRR29054901 SRR29054902 ... SRR29055188 SRR29055191
## colData names(42): Label Qubit_ng_ul ... Concentration..ng.uL. is.neg
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## rowLinks: NULL
## rowTree: NULL
## colLinks: NULL
## colTree: NULL

8.2 ANCOM II

8.2.1 Genus

set.seed(123) # Set seed for reproducibility 

# Genus level
fecal.Q2.Genus <- ANCOMBC::ancom(data=tse.fecal.Q2.Genus, assay_name="counts", tax_level="Genus",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Genus = fecal.Q2.Genus$res
res.fecal.Q2.Genus %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Genus.sig
tab.fecal.Q2.Genus.sig 

8.2.2 Family

set.seed(123) # Set seed for reproducibility 

# Family level
fecal.Q2.Family <- ANCOMBC::ancom(data=tse.fecal.Q2.Family, assay_name="counts", tax_level="Family",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Family = fecal.Q2.Family$res
res.fecal.Q2.Family %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Family.sig
tab.fecal.Q2.Family.sig 

8.2.3 Order

set.seed(123) # Set seed for reproducibility 

# Order level
fecal.Q2.Order <- ANCOMBC::ancom(data=tse.fecal.Q2.Order, assay_name="counts", tax_level="Order",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Order = fecal.Q2.Order$res
res.fecal.Q2.Order %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Order.sig
tab.fecal.Q2.Order.sig 

8.2.4 Class

set.seed(123) # Set seed for reproducibility 

# Class level
fecal.Q2.Class <- ANCOMBC::ancom(data=tse.fecal.Q2.Class, assay_name="counts", tax_level="Class",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Class = fecal.Q2.Class$res
res.fecal.Q2.Class %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Class.sig
tab.fecal.Q2.Class.sig 

8.2.5 Phylum

set.seed(123) # Set seed for reproducibility 

# Phylum level
fecal.Q2.Phylum <- ANCOMBC::ancom(data=tse.fecal.Q2.Phylum, assay_name="counts", tax_level="Phylum",
                               p_adj_method="holm", prv_cut=0, lib_cut=0, 
                               main_var="Diet", adj_formula="Week + Genotype",
                               alpha = 0.05, struc_zero=TRUE,
                               n_cl = n_cl_requested) 
## 'ancom' has been fully evolved to 'ancombc2'. 
## Explore the enhanced capabilities of our refined method!
res.fecal.Q2.Phylum = fecal.Q2.Phylum$res
res.fecal.Q2.Phylum %>% filter(detected_0.7 == TRUE) -> tab.fecal.Q2.Phylum.sig
tab.fecal.Q2.Phylum.sig 

8.2.6 Aggregate Results

tab.fecal.Q2.Genus.sig$Tax_rank <- "Genus"
tab.fecal.Q2.Family.sig$Tax_rank <- "Family"
tab.fecal.Q2.Order.sig$Tax_rank <- "Order"
tab.fecal.Q2.Class.sig$Tax_rank <- "Class"
tab.fecal.Q2.Phylum.sig$Tax_rank <- "Phylum"

rbind(tab.fecal.Q2.Genus.sig,
      tab.fecal.Q2.Family.sig,
      tab.fecal.Q2.Order.sig,
      tab.fecal.Q2.Class.sig,
      tab.fecal.Q2.Phylum.sig) -> tab.fecal.Q2

tab.fecal.Q2$Test <- "ANCOM II, fecal microbiome: no surgery, main effect Diet (structural zeroes) and week + genotype"

tab.fecal.Q2 

8.3 DESeq2

8.3.1 Genus

ds.fecal.Q2_Genus <- DESeq2::DESeqDataSet(tse.fecal.Q2.Genus, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Genus <- DESeq2::DESeq(ds.fecal.Q2_Genus)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Genus)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1238163       2007      45853         10035      10318.        5570.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    70          0           43.2
resultsNames(dds.fecal.Q2_Genus)
## [1] "Intercept"                         "Diet_High_Fat_Diet_vs_Normal_Chow"
## [3] "Week_10_vs_0"                      "Week_13_vs_0"                     
## [5] "Genotype_KO_vs_WT"                 "DietHigh_Fat_Diet.Week10"         
## [7] "DietHigh_Fat_Diet.Week13"
results(dds.fecal.Q2_Genus, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 
results(dds.fecal.Q2_Genus, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Genus.Res
fecal.Q2.Genus.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

results(dds.fecal.Q2_Genus, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

results(dds.fecal.Q2_Genus, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

results(dds.fecal.Q2_Genus, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Genus.Res) -> fecal.Q2.Genus.Res

fecal.Q2.Genus.Res$Compare_taxon <- "Genus"

fecal.Q2.Genus.Res

8.3.2 Family

ds.fecal.Q2_Family <- DESeq2::DESeqDataSet(tse.fecal.Q2.Family, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Family <- DESeq2::DESeq(ds.fecal.Q2_Family)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Family)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1571691       2593      56081        12346.      13097.        6989.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    37          0           25.3
results(dds.fecal.Q2_Family, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) -> fecal.Q2.Family.Res
fecal.Q2.Family.Res$Comparison <- c("Diet: HFD vs. NC")

results(dds.fecal.Q2_Family, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Genotype: KO vs. WT")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

results(dds.fecal.Q2_Family, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

results(dds.fecal.Q2_Family, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 10 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

results(dds.fecal.Q2_Family, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Family.Res) -> fecal.Q2.Family.Res

fecal.Q2.Family.Res$Compare_taxon <- "Family"

fecal.Q2.Family.Res

8.3.3 Order

ds.fecal.Q2_Order <- DESeq2::DESeqDataSet(tse.fecal.Q2.Order, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Order <- DESeq2::DESeq(ds.fecal.Q2_Order)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Order)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1593986       2595      57077         12560      13283.        7111.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    26          0           18.0
results(dds.fecal.Q2_Order, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 
results(dds.fecal.Q2_Order, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Order.Res
fecal.Q2.Order.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res

results(dds.fecal.Q2_Order, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res

results(dds.fecal.Q2_Order, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) 
results(dds.fecal.Q2_Order, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Order.Res) -> fecal.Q2.Order.Res

fecal.Q2.Order.Res$Compare_taxon <- "Order"

fecal.Q2.Order.Res

8.3.4 Class

ds.fecal.Q2_Class <- DESeq2::DESeqDataSet(tse.fecal.Q2.Class, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Class <- DESeq2::DESeq(ds.fecal.Q2_Class)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Class)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1594382       2595      57079         12562      13287.        7112.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    12          0           10.6
results(dds.fecal.Q2_Class, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 
results(dds.fecal.Q2_Class, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Class.Res
fecal.Q2.Class.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res

results(dds.fecal.Q2_Class, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res

results(dds.fecal.Q2_Class, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) 
results(dds.fecal.Q2_Class, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Class.Res) -> fecal.Q2.Class.Res

fecal.Q2.Class.Res$Compare_taxon <- "Class"

fecal.Q2.Class.Res

8.3.5 Phylum

ds.fecal.Q2_Phylum <- DESeq2::DESeqDataSet(tse.fecal.Q2.Phylum, ~ Diet * Week + Genotype)
## converting counts to integer mode
dds.fecal.Q2_Phylum <- DESeq2::DESeq(ds.fecal.Q2_Phylum)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq2::summary(dds.fecal.Q2_Phylum)
## $samples
## # A tibble: 1 × 6
##   total_counts min_counts max_counts median_counts mean_counts stdev_counts
##          <dbl>      <dbl>      <dbl>         <dbl>       <dbl>        <dbl>
## 1      1594436       2595      57083         12562      13287.        7112.
## 
## $features
## # A tibble: 1 × 3
##   total singletons per_sample_avg
##   <int>      <int>          <dbl>
## 1    10          0           8.79
results(dds.fecal.Q2_Phylum, tidy = TRUE, name =c("Diet_High_Fat_Diet_vs_Normal_Chow")) %>% # tidy converts to df
  filter(padj < 0.05) 
results(dds.fecal.Q2_Phylum, name =c("Genotype_KO_vs_WT"), tidy = TRUE) %>%
  filter(padj < 0.05) -> fecal.Q2.Phylum.Res
fecal.Q2.Phylum.Res$Comparison <- c("Genotype: KO vs. WT")

results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week10"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 10")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res

results(dds.fecal.Q2_Phylum, name =c("DietHigh_Fat_Diet.Week13"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("HFD, wk 13")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res

results(dds.fecal.Q2_Phylum, name =c("Week_10_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) 
results(dds.fecal.Q2_Phylum, name =c("Week_13_vs_0"), tidy = TRUE) %>%
  filter(padj < 0.05) -> tmp
tmp$Comparison <- c("Week: 13 vs. 0")
rbind(tmp, fecal.Q2.Phylum.Res) -> fecal.Q2.Phylum.Res

fecal.Q2.Phylum.Res$Compare_taxon <- "Phylum"

fecal.Q2.Phylum.Res

8.3.6 Aggregate results

rbind(fecal.Q2.Genus.Res, 
      fecal.Q2.Family.Res,
      fecal.Q2.Class.Res,
      fecal.Q2.Order.Res,
      fecal.Q2.Phylum.Res) -> fecal.Q2.deseq

fecal.Q2.deseq$Test <- "DESeq2, Fecal microbiome: no surgery, ~ Diet * Week + Genotype"

fecal.Q2.deseq

8.4 Intersection

tab.fecal.Q2 %>% filter(taxon %in% fecal.Q2.deseq$row) -> ancom.intersect.nosurgery
ancom.intersect.nosurgery
fecal.Q2.deseq %>% filter(row %in% tab.fecal.Q2$taxon) %>%
  dplyr::select(row, Comparison, Compare_taxon, Test)  %>%
  group_by(row) %>%
  mutate(n.ind = paste0(1:n()))  %>%
  pivot_wider(
    names_from = "n.ind",
    values_from = Comparison
  )  -> fecal.Q2.deseq.intersect.wide

fecal.Q2.deseq.intersect.wide

8.5 Graph

fecal.RA.Full %>% filter(Surgery == "None") -> fecal.graph.nosurgery
head(fecal.graph.nosurgery)
tail(fecal.graph.nosurgery)
for (i in 1:length(ancom.intersect.nosurgery$taxon)) {
  fecal.graph.nosurgery %>% 
  filter(OTU == ancom.intersect.nosurgery$taxon[i]) %>%  
  ggplot(aes(x = as.character(Week), y = Abundance, color = Diet)) + 
  geom_boxplot() + geom_point() +  
  labs(y="Relative Abundance", title = str_glue(ancom.intersect.nosurgery$taxon[i], " Relative Abundance, fecal"),
       x="Weel") + 
  facet_wrap(~Diet) + 
  theme_bw(base_size = 14) + 
  theme(legend.position = "top") -> p
  print(p)
}

9 Export csv files

fecal.Q1.deseq %>% filter(row %in% ancom.intersect.hfd$taxon) -> fecal.Q1.deseq.intersect
fecal.Q2.deseq %>% filter(row %in% ancom.intersect.nosurgery$taxon) -> fecal.Q2.deseq.intersect
write_csv(fecal.RA.Full, file.path(git.dir, "fecal_relative_abundance.csv"), append = FALSE, col_names = TRUE)
write_csv(ancom.intersect.hfd, file.path(git.dir, "fecal_ancom_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(ancom.intersect.nosurgery, file.path(git.dir, "fecal_ancom_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(fecal.Q1.deseq.intersect, file.path(git.dir, "fecal_deseq_hfd_intersect.csv"), append = FALSE, col_names = TRUE)
write_csv(fecal.Q2.deseq.intersect, file.path(git.dir, "fecal_deseq_nosurgery_intersect.csv"), append = FALSE, col_names = TRUE)

10 Reproducibility

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] stringr_1.5.0               tibble_3.2.1               
##  [3] readr_2.1.4                 tidyr_1.3.0                
##  [5] dplyr_1.1.3                 DESeq2_1.40.2              
##  [7] SummarizedExperiment_1.30.2 Biobase_2.60.0             
##  [9] MatrixGenerics_1.12.3       matrixStats_1.1.0          
## [11] GenomicRanges_1.52.1        GenomeInfoDb_1.36.4        
## [13] IRanges_2.34.1              S4Vectors_0.38.2           
## [15] BiocGenerics_0.46.0         plotly_4.10.3              
## [17] DT_0.29                     ggrepel_0.9.4              
## [19] microbiome_1.22.0           colorBlindness_0.1.9       
## [21] RColorBrewer_1.1-3          pals_1.8                   
## [23] ANCOMBC_2.2.2               microViz_0.11.0            
## [25] ggpubr_0.6.0                ggplot2_3.4.4              
## [27] rstatix_0.7.2               fs_1.6.3                   
## [29] phyloseq_1.44.0            
## 
## loaded via a namespace (and not attached):
##   [1] bitops_1.0-7                   DirichletMultinomial_1.42.0   
##   [3] httr_1.4.7                     doParallel_1.0.17             
##   [5] numDeriv_2016.8-1.1            tools_4.3.1                   
##   [7] doRNG_1.8.6                    backports_1.4.1               
##   [9] utf8_1.2.4                     R6_2.5.1                      
##  [11] vegan_2.6-4                    lazyeval_0.2.2                
##  [13] mgcv_1.9-0                     rhdf5filters_1.12.1           
##  [15] permute_0.9-7                  withr_2.5.1                   
##  [17] gridExtra_2.3                  textshaping_0.3.7             
##  [19] cli_3.6.1                      sandwich_3.1-0                
##  [21] labeling_0.4.3                 sass_0.4.7                    
##  [23] mvtnorm_1.2-3                  proxy_0.4-27                  
##  [25] systemfonts_1.0.5              yulab.utils_0.1.0             
##  [27] foreign_0.8-85                 dichromat_2.0-0.1             
##  [29] scater_1.28.0                  decontam_1.20.0               
##  [31] maps_3.4.1                     readxl_1.4.3                  
##  [33] rstudioapi_0.15.0              RSQLite_2.3.3                 
##  [35] generics_0.1.3                 gridGraphics_0.5-1            
##  [37] vroom_1.6.4                    gtools_3.9.4                  
##  [39] car_3.1-2                      Matrix_1.6-1.1                
##  [41] biomformat_1.28.0              ggbeeswarm_0.7.2              
##  [43] fansi_1.0.5                    DescTools_0.99.52             
##  [45] DECIPHER_2.28.0                abind_1.4-5                   
##  [47] lifecycle_1.0.3                multcomp_1.4-25               
##  [49] yaml_2.3.7                     carData_3.0-5                 
##  [51] rhdf5_2.44.0                   Rtsne_0.16                    
##  [53] grid_4.3.1                     blob_1.2.4                    
##  [55] crayon_1.5.2                   lattice_0.22-5                
##  [57] beachmat_2.16.0                cowplot_1.1.1                 
##  [59] mapproj_1.2.11                 pillar_1.9.0                  
##  [61] knitr_1.44                     boot_1.3-28.1                 
##  [63] gld_2.6.6                      codetools_0.2-19              
##  [65] glue_1.6.2                     data.table_1.14.8             
##  [67] MultiAssayExperiment_1.26.0    vctrs_0.6.4                   
##  [69] treeio_1.24.3                  Rdpack_2.6                    
##  [71] cellranger_1.1.0               gtable_0.3.4                  
##  [73] cachem_1.0.8                   xfun_0.40                     
##  [75] rbibutils_2.2.16               S4Arrays_1.2.1                
##  [77] survival_3.5-7                 SingleCellExperiment_1.22.0   
##  [79] iterators_1.0.14               gmp_0.7-2                     
##  [81] TH.data_1.1-2                  nlme_3.1-163                  
##  [83] bit64_4.0.5                    bslib_0.5.1                   
##  [85] irlba_2.3.5.1                  vipor_0.4.5                   
##  [87] rpart_4.1.21                   colorspace_2.1-0              
##  [89] DBI_1.1.3                      Hmisc_5.1-1                   
##  [91] nnet_7.3-19                    ade4_1.7-22                   
##  [93] Exact_3.2                      tidyselect_1.2.0              
##  [95] bit_4.0.5                      compiler_4.3.1                
##  [97] htmlTable_2.4.2                BiocNeighbors_1.18.0          
##  [99] expm_0.999-7                   DelayedArray_0.26.7           
## [101] checkmate_2.3.0                scales_1.2.1                  
## [103] digest_0.6.33                  minqa_1.2.6                   
## [105] rmarkdown_2.25                 XVector_0.40.0                
## [107] htmltools_0.5.6.1              pkgconfig_2.0.3               
## [109] base64enc_0.1-3                lme4_1.1-34                   
## [111] sparseMatrixStats_1.12.2       fastmap_1.1.1                 
## [113] rlang_1.1.1                    htmlwidgets_1.6.2             
## [115] DelayedMatrixStats_1.22.6      farver_2.1.1                  
## [117] jquerylib_0.1.4                zoo_1.8-12                    
## [119] jsonlite_1.8.7                 energy_1.7-11                 
## [121] BiocParallel_1.34.2            BiocSingular_1.16.0           
## [123] RCurl_1.98-1.13                magrittr_2.0.3                
## [125] Formula_1.2-5                  scuttle_1.10.3                
## [127] GenomeInfoDbData_1.2.10        Rhdf5lib_1.22.1               
## [129] munsell_0.5.0                  Rcpp_1.0.11                   
## [131] ape_5.7-1                      viridis_0.6.4                 
## [133] CVXR_1.0-11                    stringi_1.7.12                
## [135] rootSolve_1.8.2.4              zlibbioc_1.46.0               
## [137] MASS_7.3-60                    plyr_1.8.9                    
## [139] lmom_3.0                       Biostrings_2.68.1             
## [141] splines_4.3.1                  multtest_2.56.0               
## [143] hms_1.1.3                      locfit_1.5-9.8                
## [145] igraph_1.5.1                   ggsignif_0.6.4                
## [147] rngtools_1.5.2                 reshape2_1.4.4                
## [149] ScaledMatrix_1.8.1             evaluate_0.22                 
## [151] tzdb_0.4.0                     nloptr_2.0.3                  
## [153] foreach_1.5.2                  purrr_1.0.2                   
## [155] rsvd_1.0.5                     broom_1.0.5                   
## [157] Rmpfr_0.9-3                    e1071_1.7-13                  
## [159] tidytree_0.4.5                 ragg_1.2.6                    
## [161] viridisLite_0.4.2              class_7.3-22                  
## [163] gsl_2.1-8                      lmerTest_3.1-3                
## [165] memoise_2.0.1                  beeswarm_0.4.0                
## [167] cluster_2.1.4                  TreeSummarizedExperiment_2.8.0
## [169] mia_1.8.0